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Abstract 

We present in a unified framework new conforming and nonconforming Virtual Element Meth¬ 
ods (VEM) for general second order elliptic problems in two and three dimensions. The differential 
operator is split into its symmetric and non-symmetric parts and conditions for stability and accu¬ 
racy on their discrete counterparts are established. These conditions are shown to lead to optimal 
H 1 - and id-error estimates, confirmed by numerical experiments on a set of polygonal meshes. The 
accuracy of the numerical approximation provided by the two methods is shown to be comparable. 


1 Introduction 

The Virtual Element Method (VEM) was introduced in [5] as a generalisation of the conforming 
finite element method (FEM), offering great flexibility in utilising meshes with (almost) arbitrary 
polygonal elements. Unlike the polygonal finite element method (PFEM) [29] and other conforming 
FEM extensions based on the Generalised Finite Element [4] framework such as the CFE method [22] 
and the XFEM [20], the VEM handles meshes with general shaped elements in a manner that avoids the 
explicit evaluation of the shape functions. Indeed, the VEM only requires the knowledge of a polynomial 
subspace of the local finite element space to provide stable and accurate numerical methods. This feat 
is achieved by separating the contributions of the polynomial subspace from that of the remaining 
non-polynomial virtual subspace through the introduction of suitable projection operators which can 
be computed just using the VEM degrees of freedom. The polynomial consistency terms of the bilinear 
form, responsible for convergence properties of the method, are computed accurately. The remaining 
terms are only required to ensure the stability of the method, and hence they can be rougly estimated 
from the degrees of freedom. The VEM approach can also be viewed as a variational analogue of the 
mimetic finite difference (MFD) method; see [9] and the recent review paper [25]. As such, for its 
analysis we can take advantage of the standard tools of finite element analysis. 

An alternative approach is to completely relax the inter-element conformity requirements for the 
discrete space, so that simple (polynomial) spaces can be used. For instance, discontinuous Galerkin 
and weak Galerkin methods, whereby inter-element continuity is weakly imposed, are naturally suited 
to general meshes; see [15, 17, 28] and the references therein. 

Here we shall stop just short of that, presenting a general VEM framework which is based on 
relaying on some form of continuity, in the spirit of [18]. The framework is used to introduce a 
C'°-conforming and a nonconforming VEM. 

The original VEM in [5] is a (^-conforming method for solving the two-dimensional Poisson equa¬ 
tion and the same problem is considered in [3], where a nonconforming formulation is presented. The 
extension of these methods to general elliptic problems with variable coefficients in two and three di¬ 
mensions, is non trivial. Diffusion problems with non-constant diffusion tensors in two dimensions are 
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treated in [10], where VEMs which incorporate inter-element continuity of arbitrary degree are pre¬ 
sented. A crucial step towards the inclusion of low-order differential terms is provided in [1] with the 
extension of the original C-conforming VEM to reaction-diffusion problems with constant coefficients 
in two and three dimensions. This approach is extended to the solution of general elliptic problems 
in two dimensions in [8]. Concurrently, the VEM framework has been extended to the solution of 
plate-bending problems [14], linear elasticity problems in two and three spatial dimensions [6, 21], 
the Steklov eigenvalue problem [27], the simulation of discrete fracture networks [11], and the two- 
dimensional streamline formulation of the Stokes problem [2]. 

We present here a conforming and a nonconforming VEM for the numerical treatment of general 
linear elliptic problems with variable coefficients in two and three spatial dimensions. The accuracy 
and stability of the two methods is determined in Section 3 through a unified abstract framework, cf. 
Assumption Al. Here the partial differential operator is split into its symmetric and skew-symmetric 
parts and the VEM polynomial consistency and stability properties are established for each of these 
components separately, cf. Assumption A2. This approach is quite natural in that, for instance, it 
is clear that only the symmetric component is needed for the method’s stability. Indeed a unique 
stabilisation for all the terms that contribute to the symmetric part (the diffusion, reaction, and 
symmetric contribution of the convection term) is introduced. The stabilisation automatically adjusts 
with the relative magnitude of the (symmetric) terms. It also leads the way to the design of ad hoc 
stabilisation techniques for the pre-asymptotically stable solution of convection-dominated problems, 
such as the classical streamline diffusion method [23], although this is not considered here. 

To deal with non-constant coefficients and the lower-order terms, we take the approach of [8], 
rather than that of [10]. A crucial role in the former formulation is played by the L 2 -projector which 
maps the functions of the virtual element space and their gradients onto polynomials. In order to 
have the the L 2 -projection operator computable by using only the degrees of freedom, in Section 4 
we generalise a procedure introduced in [1], dubbed VEM enhancement, used here in the context of 
nonconforming VEM for the first time. In this way a family of virtual element spaces is defined from 
which a particularly simple choice can be made, cf. Section 7. This approach differs completely from 
that presented for a non-constant diffusion tensor in [10], which required the construction of a bespoke 
projection operator dependent on the diffusion tensor. The key advantage of removing this dependence 
is that lower order terms can be dealt with in an identical manner. Furthermore, we are now easily 
able to analyse the impact on the method of the approximation of the problem’s coefficients. Here 
it is important to stress that such approximation is only needed to compute the integrals involved 
in the polynomial consistency terms of the bilinear form. This fact is discussed in Section 6, with 
the conclusion that the stability and optimal accuracy of the method based on using polynomials of 
order up to k are unaffected by the use of a quadrature scheme to approximate the consistency terms, 
provided that this is of at least degree 2k — 2. We stress that this is exactly the same requirement of 
the finite element methods [16]. 

The new unified formulation offers some indisputable advantages. From a theoretical viewpoint, 
it permits us to analyse in a unified manner the conforming and nonconforming VEM following the 
standard analyses of finite element methods for elliptic problems. The analysis, detailed in Section 5, 
ultimately leads to optimal order H 1 - and L 2 -error estimates for both methods under the same reg¬ 
ularity assumptions on the mesh and the exact solution. By contrast, the analysis of conforming 
VEMs for the same problem in two space dimensions given in [8] is based on an inf-sup argument 
relaying on the mesh size being small enough. From a practical viewpoint, the implementation of the 
conforming and nonconforming VEM is formally the same (see Section 7). In fact, the only differ¬ 
ence is in the construction of the L 2 projection operator for the shape functions and their gradients. 
Such construction depends on the degrees of freedom, which necessarily differ for the conforming and 
nonconforming VEM. Also, as mentioned above, we foresee that the present unified framework will 
facilitate the treatment of convection-dominated diffusion problems, as well as the design of VEM for 
the Stokes system. 

The conforming and nonconforming VEMs are assessed in Section 8 solving numerically a rep¬ 
resentative convection-reaction-diffusion problem with variable coefficients in two dimensions. The 
accuracy of the numerical approximation provided by the two methods is comparable and confirms the 
optimal convergence rates in the L 2 - and H l - norm established by the theoretical analysis presented 
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in Section 5. Finally, in Section 9 we offer our final conclusions. 


2 The Continuous Problem 

Consider the boundary value problem 

—V • (k(i)Vu) + (3{x) ■ Vw + 'y(x)u = f(x) in (2.1a) 

u = 0 on dfi, (2.1b) 

where C lR rf is a polygonal domain for d = 2 and a polyhedral domain for d = 3. We assume that 
the coefficients Kij(x), f3i(x),'y(x) are in L°°(f2), and / £ L°°(Q) is the forcing function. We further 

suppose that k(:e) is a full symmetric d x d diffusivity tensor and is strongly elliptic, i.e. there exist 

k*, k* > 0, independent of v and x, such that 

n/v(x)\ 2 < v(x) ■ k(x)v(x) < k*|u(£c)| 2 , (2.2) 

for almost every x £ tt and for any v £ (7fo(fl)) d , where |-| denotes the standard Euclidean norm on 
]R d . Finally, we suppose that there exists /x 0 > 0 such that 

M*) : = 7(*) - • (3{x) > no > 0, (2.3) 

for almost every x £ Cl, and assume that V • (3 £ 

The variational form of problem (2.1) reads: find u £ such that 

(kVu, Vv) + (/3 • Vm, v) + (ju, v) = (/, v) \/v £ (2.4) 

with (•, •) denoting the L 2 inner product on f2. We split the bilinear form on the left-hand side of (2.4) 

into its symmetric and skew-symmetric parts: 

a(u, v) := (kVm, Vv) + (/xu, v ), (2.5a) 

b(u, v) := ^ [{(3 ■ Vm, v) - (u, f3 ■ Vw)], (2.5b) 

and consider discretising the problem written in the equivalent form: find u £ Hq(Q) such that 

A(u,v) := a(u,v) + b{u,v) = (f,v) \/v £ Hq(CI). (2-6) 

Rewriting the variational form in this way would not be necessary for the classical finite element 
method, but it turns out to be a useful step for the Virtual Element Method in view of ensuring that 
the discrete framework preserves the information about the symmetric and skew-symmetric parts of 
the bilinear form. 

It is simple to check that the bilinear form A is coercive and bounded, and the variational problem 
therefore possesses a unique solution by the Lax-Milgram lemma. 


3 The Virtual Element Framework 

We assume that a Virtual Element Method (VEM) consists of the following fundamental ingredients: 

Assumption Al. For any fixed h > 0 and k £ N, we have: 

• A finite decomposition (mesh) {7/,} of the domain Q into non-overlapping simple polygonal/ 
polyhedral elements with maximum size h. The adjective simple refers to the fact that the 
boundary of each element in the decomposition must be non-intersecting. Further, the boundary 
of any element E £ 77 is made of a uniformly bounded number of interfaces (edges/faces) which 
are either part of the boundary of f1 or shared with another element in the decomposition. 
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• A finite dimensional function space Vj t C H 1 (Th) where 

H\T h ) := {n G L 2 (H) : v\ E G G T4 (3.1) 

(not necessarily a subspace of H ( \ (Q)) to be used as a trial and test space, on which the following 
Poincare-Friedrichs inequality for piecewise H 1 functions holds: 

Iloilo,a — ^pf|^/ i| i^h =: ^2 Iloilo,(3-2) 

EGTh 

hence | - |i,/» is a norm on V/,. 

Further, for each element G Th, the space V E := Vh\E must contain the space Vk{E) of 
polynomials of degree k on E ; 

• A bilinear form Ah : Vh x Vh —t IR, which may be split over the elements in the mesh Th as 

A h (u h ,v h )= ^ A E (u h ,v h ), 

E&T h 

for any Uh,Vh G Vh, where A E is a bilinear form over the space V E . 

• An element //, G approximating the forcing term. 

Remark 1. The definition of simple polygons and simple polyhedra is general enough to include, 
for instance, elements with consecutive co-planar edges/faces, such as those typical of locally refined 
meshes with hanging nodes and non-convex elements. Later on, in Assumption (AS) in Section f.f, we 
shall introduce some standard mesh regularity assumptions which are required for the approximation 
properties of the virtual element spaces. 

In view of the following analysis, it is useful to extend the definition of the continuous bilinear 
form A, as well as of its parts a and b , to the whole of H l (Th) as a sum of elemental contributions. 
Henceforth, we shall use 

A(u,v) := ^2 A E (u,v)= ^2 a E (u,v) + b E (u,v) Vu,v G H l (Th). 

EClTh EeTh 

We also place a few more restrictions on the nature of the bilinear forms A E . As with the continuous 
bilinear form, we write A E as the sum of a symmetric and a skew-symmetric part: 

A E (u h ,v h ) = a E (u h ,v h ) + b E (u h ,v h ), 

and require that they satisfy the following properties: 

Assumption A2. The bilinear forms a E and b E are assumed to satisfy the properties of polynomial 
consistency and stability, defined as 

• Polynomial consistency: If either uu G Vk{E) or Vh G Vk{E), the symmetric and skew-symmetric 
parts of the local virtual element bilinear form must satisfy 

a E {u h ,v h )= [ KU 0 k _ 1 (Tuh)-U 0 k _ 1 {yv h )dx+ f pU° k u h Il 0 k Vhdx, 

J E J E 

b E iuh,v h ) = \J e @' [ntitv^nK -nguX-itVwfc)] dx, 

where the operator n° : L 2 (E ) —> Ve(E) for £ < k denotes the L 2 (£')-orthogonal projection onto 
the polynomial space Ve(E), and is defined for any function v G L 2 (E) as the unique element 
n°u of Vi(E) such that 

(n ° e v,p) E = ( v,p) E Vp G Vi{E). 
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• Stability: There exist positive constants a*, a*, and (3* independent of h and the mesh element 
E such that, for all Vh,Wh G V E , the symmetric part satisfies 

a*a E (v h ,v h ) < a E (v h ,v h ) < a*a E (v h ,v h ), 

and the skew-symmetric part satisfies 

b E {v h ,v h ) = 0 and b E {v h ,w h ) < ^*\\v h \\ 1E \\w h \\ 1E . 


Precise methods of choosing the spaces and bilinear forms of the method are, of course, the focus 
of most of the remainder of this paper. However, the simple properties presented above are enough 
to prove two crucial facts about the behaviour of such a Virtual Element Method: firstly that any 
such method possesses a unique solution and secondly an abstract Strang-type convergence result, 
which will be used later on to derive optimal order error bounds in the H 1 norm. These results are 
encapsulated in the following theorems. 

Theorem 1 (Existence and uniqueness of a virtual element solution). Under Assumptions Al and 
A2, the problem: find Uh G Vh such that 


ah(u h ,v h ) + b h (u h ,v h ) = (. fh,v h ) Mv h G V h , 


(3.3) 


possesses a unique solution. Here (•, •) denotes the duality pairing between Vh and its dual Vf. 

Proof. The stability condition on a E ensures that this bilinear form inherits the coercivity of its 
counterpart a E (•,•), hence 

a E (v h ,v h ) > a*(K*H v ^Ho,B + Ho\\v h \\l E ). (3.4) 

Summing up the contribution of all elements and using (3.2) we deduce the coercivity bound 

a h (v h ,v h )> 1 + ^ KHi.fc + («* + tt))IKIIn) > i + ^ ( 3 - 5 ) 

From inequality (3.4) and the symmetry and bilinearity of a E , it follows that a E is an inner product 
on V E . Hence, we can apply the Cauchy-Schwartz inequality and use the right stability inequality of 
Assumption A2 to prove the continuity of a E : 


a h(u h ,v h ) < {a%(u h ,u h ))z(a%(v h ,v h ))i < a*{a E (u h ,u h ))^{a E {v h ,v h ))^ 


<a* max{K*,||/i|| }||u h | 


l,E\\ Vh \\l,E’ 


and the continuity of ah easily follows. 

The stability property for b E means that this term does not feature in the coercivity analysis and 
imposes continuity with constant f3*. Hence we may conclude that the problem (3.3) admits a unique 
solution by the Lax-Milgram lemma. □ 


Theorem 2 (Abstract a priori error bound). Under Assumptions Al and A2, there exists a constant 
C > 0 depending only on the coercivity and the continuity constants such that 


||u 


u h\\i,h ~U 


I inf 

I Vh&Vh 


||u 


v h\\i, h 


+ sup 

Wh&Vh 

ivhAo 


l( fh,w h ) - (f,w h )I 


+ inf 

p€Vk(Th) 


||w 


p\\i, h + sup 

EdT h w h£Vff 
whAO 


A E (p,w h ) - A E (p, w h ) 


+ sup 

WhEVh 

whAo 


I A{u,w h ) - (f,w h ) 1 

\\ w h\\l,h 


(3.6) 


5 



where 


V k (T h ) := {p € L 2 (n) : p\ B G P fc (B) WE G 77,}• 

The last term in the right-hand side of the above error estimate measures the nonconformity error, i. e. 
it is non-zero only when Vh is a non-conforming virtual space. 

Proof. Let v k be an arbitrary element of 14 and let w k = «/,. — v k G Vj,. Then, the coercivity of A k 
implies that 


Oi\\u h ~ V h \\\ h < A h (u h - V h ,W h ) = (f h ,W h ) + Y [ A h (P - v h, w h) - A h (P’ W h)) 


EGTh 


= ( fh, W h ) + Y [ A h (V - Vh, Wh) - aE {p, W h ) + ( aE {p , w h ) - Af ip, to h ))] 
EaT h 

= ifh,W h ) - A(u,w h ) + Y [ A h(p - Vh,w h ) + A E (u-p,w h ) 

EeTh 

+( AE (p,w h ) - A%(p,w h ))] , 

for any p G Vk{Th). We express the potential nonconformity of the virtual element space Vh as 

A(u,w h ) = (f,w h ) + ( A(u,w h ) - (f,w h )), 


so that 


a\\ u h - v h \\\ h < [(fh, w h ) - (/, w h )] + [(/, w h ) - A[;u , w h )\ 

+ [ A h(p- v h,Wh) + A E (u-p,w h ) + (A E (p,w h )-Af l (p,w h ))]. 
E&T h 


Hence, for all Wh G Vh \ {0} and p G V k (J 7), applying the continuity of the bilinear forms, multiplying 
and dividing by ||ti;/ l || 1 h , and using the triangle inequality we find that 


a 


Wh-v h \\ l h < 


\(fh,w h ) - (f,w h )I |(/, w h ) - A(u,w h )\ 


\ w h\\ lth 


\\ w h\\ lth 


+ \\ U - V h\\ 1 , h + 2 \\ U -P\\i,h+ 

EaTh 


AE (p,w h ) - A E (p,w h ) 

II w ?i|| 1 ,E 


The result now follows by the triangle inequality. 

4 The Virtual Element Spaces 

We introduce two types of spaces in two- and three-dimensions implementing the framework of Sec¬ 
tion 3: a conforming and a nonconforming virtual element space. With deliberate ambiguity we refer 
to both spaces as 14 to emphasise the fact that the method is otherwise the same in either case. In 
all cases the local virtual element space V E must contain the space Vk(E) of polynomials of degree up 
to k on E, cf. Assumption Al. The complement of Vk(E) in V E is made up of functions which are 
deemed expensive to evaluate, although they may be described through a set of (known) degrees of 
freedom. Central to the virtual element methodology is the idea of computability , defined as follows: 

Definition 1. A term will be called computable if it may be evaluated using just the degrees of freedom 
of a function and the polynomial component of the virtual element space. 

We require that the polynomial consistency forms of Assumption A2 are computable. Hence the 
spaces must be constructed in such a way that the projections 

n^Vi^ and (4.1) 
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are computable for any function v k £ V/f. While the first of these is computable for functions in the 
original conforming (for d = 2) and nonconforming virtual element spaces presented in [5, 3], the second 
is not. However, a modification to the conforming space was introduced in [f] which renders n^Vh 
computable without changing the degrees of freedom used to describe the space. For the construction 
of the spaces which follow, we consider a generalisation of the process presented in [1], 

We first introduce an appropriately scaled basis for V k (E ), k £ N. Denote by A i}(E), feM, the 
set of scaled monomials 


M}{E) 


X — Xe 

He 


, |s| =(- 


where a is a multi-index with |s| := sj + s 2 and x s := x^x^ 2 . Further, we define Ai k (E) := 
\Ji< k M*(E) {m a }a=i, a basis of V k {E), where N d , k := dim(V k <JR d ))- 

Below, s denotes a d — 1 dimensional mesh interface (either an edge when d = 2 or a face when 
d = 3) of the mesh element E, and the set of all mesh interfaces in 77, will be denoted by <S>/,. This set is 
divided into the set of boundary edges <Sj( dly := {s £ S k : s C 9H} and internal edges S ^ nt := Sh\S^ diy . 
Bases for polynomial spaces defined on a interface s can be similarly constructed; the same notation 
will be used. We denote by ve the number of interfaces s £ dE. 


4.1 The Local Spaces 

The construction of the local nonconforming virtual element space is formally the same for d = 2 or 
3, while the construction of the conforming space is hierarchical in the space dimension. Because of 
this, while we simultaneously provide here a definition of the nonconforming space for d = 2 or 3, we 
initially only consider the conforming space for d = 2, postponing the construction for d = 3 until the 
end of the section. 

We first introduce the two sets of degrees of freedom used in [5, 3] to describe the original local 
virtual element spaces. The final spaces presented later in the section will be described using exactly 
the same degrees of freedom, so we record them here. 

Definition 2. The degrees of freedom for the local conforming and nonconforming spaces are 
(a) for the conforming space 

• the value of vh at each vertex of E; 

• for k > 1, the moments of Vh of up to order k — 2 on each mesh interface s C dE 

v h m a ds Vm a £ At fc _ 2 (s); 

for the nonconforming space, the moments ofvh of up to order k — 1 on each mesh interface s C dE 

v h m a ds V?n a £ M k -i(s ); 

(■ b) for k > 1, the moments of Vh of up to order k — 2 inside the element E 

tL / v h m a dx \/m a £ M k - 2 {E). 

\E\ Je 

A counting argument shows that the cardinality of the above sets of degrees of freedom is ue = 
ve + t'EA r i,t -2 + -/V 2) fc_ 2 for the conforming case and n ,e = VENd-i tk -i + for the nonconforming 
case. The degrees of freedom for a hexagonal element are represented in Figure 1 and Figure 2 for the 
nonconforming and conforming case, respectively. The nonconforming degrees of freedom for a cubic 
element are shown in Figure 3. 
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k = 1 


k = 2 


k = 3 


k = 4 


Figure 1: Degrees of freedom of the non-conforming VEM for a hexagonal mesh element for k = 
1, 2, 3,4; edge moments are marked by a circle; internal moments are marked by a square. 






Figure 2: Degrees of freedom of the conforming VEM for a hexagonal mesh element for k = 1,2,3,4; 
vertex values and edge moments are marked by a circle; internal moments are marked by a square. 


We first consider a conforming and a nonconforming superspace in which both terms in (4.1) are 
computable. Again with deliberate ambiguity, we refer to the enlarged space in each case as . For 
the enlarged nonconforming space on the element E, we define 


O 

W* -=^h& H\E) : Av h £ V k {E) and £ P fc _i(s) Vs C dE 


For the conforming space, we first introduce the boundary space 


B k (dE) := {u £ C°(dE) : w| s £ V k (s) for each interface s of dE } , 
and define as 


Wf := {v h £ H 1 ^) : Av h £ V k (E) and v h \ dE £ B k (dE )} . 


It is clear that, in either case, we have V k {E) C W/f. Further, the space W^f may be described using 
the combination of the degrees of freedom of Definition 2 and the extra degrees of freedom. 


Definition 3. The extra degrees of freedom are taken as the moments of v k of order k and 
inside the element E 


1 

W\ 


Vhm a dx 


Vm a £ M%(E) U Ml-^E). 


k - 


1 


The proof that the combined set of degrees of freedom is unisolvent for the conforming space is given 
in [1], and for the nonconforming space follows exactly as the original unisolvence proof in [3]. It is 
based on the observation that each space defines its elements as those functions which solve a particular 
class of Poisson problem, with piecewise polynomial Dirichlet and Neumann boundary conditions in 
the conforming and nonconforming case, respectively, specified by the degrees of freedom. Similarly, 
we can easily prove the crucial fact that any p k £ V k (E) is uniquely determined by the original degrees 
of freedom of Definition 2, cf. [5, 3]. Indeed, if p k £ V k (E) and all the original degrees of freedom of 
p k are zero then, 

f dp 

(Vp fc , Vp fe ) E = (-Ap k p k ) E + / — Pfcds = 0. 

j dE o n 
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The first term is zero as A p k = 0 if k = 1 and as A p k £ Vk-2 (E) if k > 1 and hence the first term is 
a linear combination of the (zero) internal degrees of freedom of pk in this case. The second term is 
also zero because it is always a linear combination of the boundary degrees of freedom of p k . Hence, 
Pk = constant in E. The fact that p k = 0 now follows from the hypothesis that pk is zero at any 
vertex in the conforming case, and f s Pk ds = 0 for any s £ dE in the nonconforming case. 

Using collectively the degrees of freedom of Definitions 2 and 3, it is possible to compute both of 
the projections in (4.1) for any Vh £ Wjf. Calculating nJJ'Ch requires solving the variational problem: 
find £ Vk(E) such that 

(n° k v h ,m a ) E = {v h ,m a ) E \/m. a £M k (E), (4.2) 

which is computable as the quantities on the right-hand side are the internal degrees of freedom of Vh- 
Similarly, computing requires finding the polynomial W k _ 1 '\7vh £ (Vk{E)) such that 

{Hk-i'^v h ,m a ) E = (Vv h ,m a ) E Vm a £ (Mk-i(E)) 2 . (4.3) 

This is possible because the right-hand side is a sum of boundary and internal degrees of freedom of 
Vh'- 

(Vv h ,m a ) E = n • m a v h ds - (v h , V • m a ) E . 

J dE 

Note that this time only the original degrees of freedom of Definition 2 are required. 

However, in each case this enlarged space requires an extra card(Al^(£ ; )) + card(Al))_ 1 (U)) degrees 
of freedom (namely the extra internal moments in Definition 3) compared with the original spaces 
introduced in [5] and [3], which are described by the degrees of freedom of Definition 2 only. To 
reduce the number of degrees of freedom, we adopt a generalisation of the procedure introduced in [1], 
producing a family of different subspaces of W)f spanned by the original sets of degrees of freedom 
of Definition 2, yet in which we can still compute the required projections in (4.1). The procedure 
consists of the following three steps. 

1) We introduce an equivalence relation ~ on W)f, defining Vh ~ Wh if all of the original degrees of 
freedom of Vh and Wh are equal, and consider the quotient space /~ which, by construction, 
is spanned by the original degrees of freedom of Definition 2. 

2) Since Vk{E) C W® and any p k £ Vk(E) is uniquely determined by the original degrees of 

freedom, we may conclude that any equivalence class [i>h] contains at most one polynomial. 
Then, we may unambiguously associate Vk{E) with the resulting ‘polynomial’ subspace of W)f /^. 
Hence, we can introduce any projection operator nj : —> Vk{E) C Wjf / ^ which associates 

a polynomial to each equivalence class [i>/,]. 

3) The local virtual space V k is then defined by selecting a specific representative from each equiv¬ 
alence class in For each [vh\ € W k /A we take the function Wh G [■%] such that the extra 

degrees of freedom (in Definition 3) of wh are equal to those of njlfuh]. Note in particular that 
Pk(E)cV h E . 

Remark 2. This is a generalisation of the idea introduced in [1], where only the H 1 (E)-orthogonal 
projection of Vh into Vk{E) was considered for H* . The space resulting from this choice is well defined 
because the H 1 ^)-orthogonal projection of Vh is computable using just the original degrees of freedom. 
However, the freedom in choosing n* is something we wish to exploit to produce a more computationally 
efficient method, particularly when d = 3. We explore more possible choices of n* in Section 7, although 
for now we leave this choice open. 

In more concrete terms, given a projector n£, we define the local virtual element spaces to be 

V h E := [v h £ Wf : {v h - H * k v h ,p) E = 0 Vp € M* k {E) U M * k _,(£)}, (4.4) 

1 Alternatively, one may think of this as associating a polynomial to each combination of the original degrees of 
freedom in Definition 2, independent of the extra degrees of freedom in Definition 3. 
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Figure 3: Degrees of freedom of the non-conforming VEM for a cubic mesh element for k = 1,2,3,4; 
face moments are marked by a hexagon; internal moments are marked by a square. Only the internal 
degrees of freedom and those of the visible faces are marked; the numeric labels indicate the number 
of degrees of freedom when they are more than 1. 


where W/f' denotes either the enlarged conforming or nonconforming space. Clearly, we can use the 
original degrees of freedom of Definition 2 to describe . 

Computing II ^Vh for each G is now possible, since the terms on the right-hand since of (4.2) 

are either degrees of freedom of v k or moments of II * k v h- 

4.2 The Global Spaces 

The global virtual element space in each case is constructed as a subspace of an infinite dimensional 
space V, defined differently for the conforming and nonconforming methods. For the conforming 
method, we simply take V := Uq(H). For the nonconforming method, we introduce the subspace 
Hl’ nc (T h ) of the nonconforming broken Sobolev space H 1 (7)>i) defined in (3.1), by imposing certain 
weak inter-element continuity requirements such that 

V := Hl' nc (T h ) = jw G H\T h ) : J [u] • n s q ds = 0 Vq € V k . x (s), Vs G S h 

The jump operator [■] across a mesh interface s G Sh is defined as follows for v G FI 1 (7^). If s G 5)"*, 
then there exist E + and E~ such that s C dE + n dE~. Denote by the trace of v\ ± on s from 
within £, and by nf the unit outward normal on s from E ± . Then, [v] := v + nf +v~nj. If, on the 
other hand, s G <5 ; *( dly , then [u] := vn 3 , with v representing the trace of v from within the element E 
having s as an interface and n s is the unit outward normal on s from E. 

It may be seen that the broken Sobolev norm | • |i^ is a norm on H^’ nc (Th) and hence the same 
will be true for any virtual element subspace Vh, as required by Assumption Al, cf. [12]. 

Finally, the global space is constructed in either case from the local spaces presented above as 

V h := {v h GV :v h \ E Gl)f VEeT h }. 

As global degrees of freedom we take the equivalents of those in Definition 2, namely, for each 
function Vh G Vjf, 

(a) for the conforming space 

• the value of Vh at each internal vertex of Tk’, 

• for k > 1, the moments of Vh of up to order k — 2 on each mesh interface s G 5™* 

v h m a ds Vm a G A4 fc _ 2 (s); 

for the nonconforming space, the moments of Vh of up to order fc — 1 on each mesh interface s G <S], nt 

v h m a ds \/m a G Mk-i(s)-, 
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Figure 4: Degrees of freedom of the conforming VEM for a cubic mesh element for k = 1, 2, 3,4; vertex 
values and edge moments are marked by a circle; face moments are marked by an hexagon; internal 
moments are marked by a square. Only the internal degrees of freedom and those of the visible faces 
and edges are marked; the numeric labels indicate the number of degrees of freedom when they are 
more than 1. 


(6) for k > 1, the moments of Vh of up to order k — 2 inside each element E G Th 


1 

W\ 


Vhm a dx 


Vto q G M k -2(E). 


The local degrees of freedom corresponding to boundary vertices and edges s G iS ddry are fixed as zero 
in accordance with the definition of the ambient spaces. The unisolvency of these degrees of freedom 
follows from the definition of the relevant ambient space in each case and the unisolvency of the local 
degrees of freedom. 


4.3 The Conforming Space for d = 3 

The construction of the local conforming virtual element space for d = 3 is based recursively on the 
space just detailed for d = 2. Define V® E as the 2-dimensional conforming virtual element space of 
order k constructed over the polygonal interfaces making up dE. Then, we define the space W E in 
this case to be 

Wjf := {v h G H\E) : v h \ dE G Vf E and Av h G V k (E)}. 

The degrees of freedom that we take for each function Vh G V E are then 

(а) the degrees of freedom of Vjf E ; 

(б) for k > 1, the moments of Vh of up to order k — 2 inside the element E 


1 

W\ 


Vhm a dx 


Mm a G M k -2(E). 


plus the extra degrees of freedom of Definition 3. 

This allows us to construct the space V E from Wjf in exactly the same manner detailed above. 
As before, the space V E is , spanned by just the first sets of degrees of freedom given in (a) and (b) 
above. The proof that these degrees of freedom are unisolvent is again given in [1]. Also, it is clear 
that the dimension of the local space for d = 3 is n E = v E + v' E N\,k -2 + v E N 2 ,k -2 + N ^^-2 where 
v E and u' E denote, respectively, the number of vertices and edges of E. The degrees of freedom for a 
cubic element are shown in Figure 4. 

Computing is just the same as for d = 2, since the terms on the right hand since of (4.2) are 
either degrees of freedom of Vh or moments of n^i^. To compute n^_ 1 Vt'*,, we must compute the face 
terms in (4.3). Using the L 2 (s)-orthogonal projection on the face s, these may be rewritten as 

n ■ rriall^Vh ds Vs C dE. 
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The face projection n°’ s Uh, is computable using the degrees of freedom of Vh on the face s since 
Vh\s G Vfi, and consequently this term is also computable. This means that the degrees of freedom in 
this space allow us to compute both of the required terms in (4.1). 

Finally, the global space and the set of global degrees of freedom for d = 3 are constructed from 
the local ones in the obvious way, completely analogously to the case for d = 2. 

4.4 Approximation Properties 

Both the conforming and nonconforming spaces presented above satisfy optimal approximation results 
for the approximation of sufficiently smooth functions. Since these results will be used throughout the 
remainder of the paper, we collect them together here. They rely on the following assumption on the 
regularity of the mesh Th- 

Assumption A3. (Mesh regularity). We assume the existence of a constant p > 0 such that 

• for every element E of Th and every interface s of E, h 3 > ph-E 

• every element E of Th is star-shaped with respect to a ball of radius phE 

• for d = 3, every interface s € Sh is star-shaped with respect to a ball of radius ph 3 . 

Theorem 3 (Approximation using polynomials). Suppose that Assumption A3 is satisfied. Let E £ Th 
and let 11° : L 2 (E) —> Ve(E), for £ > 0, denote the L 2 (E)~ orthogonal projection onto the polynomial 
space Vt(E). Then, for any w £ H m (E), with 1 < m < l + 1, it holds 

Ik - n °Hlo,.E + - n?Hi ,E - Ch E\ W \m,E- 

Let s be an interface shared by E + ,E~ £ Th and let II°’ S : L 2 (s) —» Vi(s), for l > 0, denote the 
L 2 (s) -orthogonal projector onto the polynomial space Vi(s). Then, for every w £ H m {E + U E~), with 
1 < m < £ + 1, it holds 


n 0,s , 7 ttO,S 

£ W o s + IT'S w — 11^ W 


, < Ch\ 

1,S — 


i—l/2| 


m,E+UE~ ’ 


In both instances, the positive constant C depends only on the polynomial degree £ and the mesh 
regularity. 

This theorem may be proven using the theory in [13] for star-shaped domains and its extension to 
more general shaped elements presented in e.g. [19]. We also have the following result regarding the 
approximation of sufficiently smooth functions by those of the virtual element space, which may be 
proven as in [27]. 

Theorem 4 (Approximation using virtual element functions). Suppose that Assumption A3 is satisfied 
and let Vh denote either the conforming or nonconforming virtual element space. Let m be a positive 
integer such that 1 < m < k + 1. Then, for any w £ there exists an element w\ £ Vh such that 


Ik - w i\\o + h \ w - wili < Ch m \w\ m 


where C is a positive constant which depends only on the polynomial degree k and the mesh regularity. 


5 Error Analysis 

We are now in a position to prove an optimal order error bound in the H 1 - and L 2 -norm for any of 
the VEM introduced above, starting with an estimate of the nonconformity error introduced by using 
the nonconforming virtual element space. 
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Lemma 1. Suppose that Assumptions A1-A3 are satisfied and let u £ H m+1 (D) for some positive 
integer m > 1 be the solution to (2.4). Define r = min (k,m) and suppose that the coefficients n,(3,n £ 
W r+1,00 (£l). Then, there exists a positive constant C independent of h and u such that 


sup 

Wh&Vh 


I A(u,w h ) - {f,w h )\ 

IMIi.h 


< Ch r \\u\\ 


r+1 ’ 


where Vh is the nonconforming virtual element space described in Section f. 


Proof. We apply the definition of the jump operator [•] and the Green’s identity under the assumption 
that m > 1, which implies that u £ and use the fact that 14 C H k nc (Th) to obtain 



- ^uf3) ■ [wj ds 

- \u(3)) • ([wh] - n°’ s lw h j) ds 


Using the Cauchy-Schwartz inequality and then applying the approximation estimates of Theorem 3 
to bound each of the resulting terms, we obtain, cf. [3] or [18], 

\A(u,Wh) ~ {fi w h)\ 4 Ch IMIr-in^+Ui;- \ w h\i t E+UE~ t 

where for each side s the symbols E + and E~ denote the two elements sharing that side, and conse¬ 
quently the lemma holds. □ 


Theorem 5 ( H 1 error bound). Suppose that Assumptions A1-A3 are satisfied. Let k > 1 be a positive 
integer and let u £ H m+1 (Q) be the true solution to problem (2.4) for some positive integer m. Define 
r = min (k,m) and suppose that the coefficients £ lU r+1,00 (U) satisfy (2.2) and (2.3). Let 

ifh,v h ) ■= YjE£T h (fh, v h)E, with f h \ E ■= n ^ ax(fc _ 2 0) /|_E. Denote by u h £ V h the corresponding 
virtual element solution to problem (3.3) where Vh is either the conforming or the nonconforming 
virtual element space presented in Section f. Then, there exists a constant C independent of h and u 
such that 

ll« - u h \\i < Ch r {\\u\\ r+1 + ||/|| T ._ 1 ) 

Proof. We prove the theorem by separately bounding the terms of the Strang-type abstract convergence 
result of Theorem 2. The first term on the right-hand side of (3.6), i.e. inf^gy^ ||tt — h , is easily 
bounded by introducing any interpolant u\ £ Vh of u as in Theorem 4. For the second term, first we 
suppose that k > 2. In this case we may apply the definition of the L 2 -projection of / to find that 


sup 

Wh&V h 


I (fh,Wh) - (/, W h )\ 

IKIli.fc. 


sup 

WhtVh. 



( n fc- 2 /~ f)wh da: 


IMIi.h 


E EaT h [ - f)(.w h -U° 0 w h )dx 

= sup --- 77 - 77 - 

w h ev h \\ w h\\i t h 

YjE^Th ll n fc—2/ _ /llo.sll^ “ n O w ft|lo,E 

< sup - 77 - 77 - 

w h ev h h 

< C^ll/llr-r 


having used the bounds of Theorem 3. A similar argument applies with fh = ITq/ when k = 1. 
Turning now to the infimum over the polynomial subspace, we first observe that 


inf 

p&Vk (Th) 


\\ u ~p\\i,h+ sup 

EeTh WheV h 


\A E {p,w h ) - A%(p,w h ) 


< Hu-ngull^ 


E sup 

BeT h w ^ v h 


\\ W h\W,E 

\A E {U° k u,w h ) ~ A E (U° k u,w h )\ 

IK^II 1,E 
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To treat the second term on the right-hand side, we recall the splitting of the bilinear forms into 
their symmetric and skew-symmetric parts, and bound each part separately. Here we only detail 
the bounding of the difference between the symmetric parts since the skew-symmetric parts can be 
treated analogously. From the definition of a E and the polynomial consistency property of a E (cf., 
Assumption (A2)) it follows that 


infeW^-n^^da; 


IE 


(!- n fc)( 7 n k u)w h dx 

0 „,\ || 


\a E {U° k u,w h ) - a%(n° k u,w h )\ < [ K\7U° k u- (I—n°_ 1 )Vw ft ,dx 

Je 

= f (i-ng^jj^vnguj-vwhd® 

Je 

< (||(i- n L 1 )('«vn°u)|| 0£; + : [j(7-n°)(7n , fe ^ ll0 ^, 

and now the results of Theorem 3 and the regularity assumption on k and /i implies 

\a E {U° k u, w h ) - af (n ° k u, w h ) 


sup 

w h ev h E 


\ w h\\i e 


< Ch E \\u\\ r+1E . 


A similar bound holds for \b E (H k u,Wh) — b E (n°u, Wh) | due to the regularity assumption on (3. Com¬ 
bining these bounds we therefore obtain 

\A E {U° k u,w h ) - A E (Il 0 k u,w h )\ 


sup 

w h ev E 


\ w h\\i,E 


< Ch r E \\u\\ r+1E . 


The result then follows by observing that Lemma 1 provides an optimal order bound for the remaining 
term relevant to the nonconforming case only. □ 

Theorem 6 (A 2 error bound). Suppose that Assumptions A1-A3 are satisfied, and further assume 
that the domain fi is convex. Let k > 1 be a positive integer and let u € be the solution to 

the problem (2.4) for some positive integer m. Define r = min(fc,m) and suppose that the coefficients 
n,f3,p& H rr+ 1 ’°°(H) satisfy (2.2) and (2.3). Let ( fh,v h ) := Y^EeT h (fh> v h)E, with f h \ E := Hg_ 1 /| jE . 
Denote by Uh G Vj, the corresponding virtual element solution to problem (3.3) where V^ is either the 
conforming or the nonconforming virtual element space presented in Section f. Then, there exists a 
constant C independent of h and u such that 


\u-u h \\o < Ch r+1 1 


llr+1' 


Proof. Let if € H' 2 (TL) n Hq(LI) be the solution to the dual problem 

—V • (nVif) — j3 ■ VV’ + (7 — V • /3)ip = u — Uh- 
Then, due to the convexity of Q, ip satisfies the regularity bound 

IHI2 < c \\ u ~ u h\\o, 

and consequently for any interpolant ipi as in Theorem 4, we have 

IIV’ - V’iII i. h < Ch\m 2 < Ch\\u - u h || 0 . 

Multiplying (5.1) by u = Uh and integrating, we find that 

||w - u h \\l = (u- u h , -V • (nVip) - f3 ■ Vip + (7 - V • f3)ip) 


(5.1) 


= A(u - Uh,ip) + y: f (\ ■ I u - u hi ds, 

ses h Js ^ ' 
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and the edge-wise term may be bounded by arguing as in Lemma 1, to find that 




• [« 


Uhj ds < C7 i||V|| 2 ||« 


U h\\l,h 


< Ch r+1 \\u\\ r+1 \\u - u h \\ 0 . 


To bound the other term, we add and subtract appropriately, 


A{u - u h , V) = A(u - u h , V - Vi) + A(u - u h , Vi) 

= A(u - u h , ip-ipi) + ( A(u , Vi) - (/, Vi)) + 

+ ( A h (u h , Vi) - A(u h , Vi)) + ((/, Vi) - (fh, Vi)), 


obtaining terms which relate to those of the original abstract error bound. We label these as Ti, T 2 , 
X 3 and Tf respectively, and bound them separately. 

We first observe that Tf may be bounded using the continuity of the variational form and the 
.fV-norm error bound of Theorem 5 as 


Ti := A(u Uhi V Vi) < C\\u - u h \\ l h \\tp - tpi\\ l h 
< Ch r+ 1 \\u\\ r+1 \\u - u h \\ 0 . 

The term T 2 measures the nonconformity of the method and may be bounded using Lemma 1 as 

\T 2 \ := \A(u,tpi)~ (/, Vi)I < C7i r ||u|| r+1 ||V ~ Villi,/* 

<Ch r+ 1 \\u\\ r+ 1 \\u-u h \\ 0 . 

Using the definition of the L 2 -projection, we can rewrite Tf as 

Tf : = (/, Vi) - (fh, Vi) = £ (/ - n ° k _J, Vi)e ='£(/- K-J, V’i - <*Pi)e 

E(z.l~h E£l~h 

^ E l|/- n °-l/|lo i B : ||^I- n oV’l||o >£ 

E£T h 

<CY h r \\f\\r, E HH\l,E<Ch r+1 E 11 / 11 ^ 11 ^ 112 ,E 

E(zil~h E^fTh 

<Ch r+ 1 \\f\\ r \\u-u h \\ 0 . 

Finally we turn to the inconsistency term T 3 , namely 

T 3 : = A h (u h , Vi) - A(u h , Vi) = E A h ( u h , Vi) - A E (u h , Vi) 

E&Th 

= E ( A h (' u h ~ n ° u > V’i - n?v) - A E (u h - n£w, Vi - n°v)) + 

EeTh 

+ (Ah (n ° k u, Vi) - A E (U° k u, Vi)) + (A E (■ u h , II? V) - A E (u h , II? V)) 

The first difference can then easily be bounded using the fact that both the variational form and the 
VEM bilinear form are continuous in the H 1 norm, so 

A E \u h - n°w, Vi - n?v) - A E (u h - u° k u, v - n?v) < c\\u h - n°u|| 1)E ||Vi - n?v|| 0)E 

< Ch r+ 1 \\u\\ r+hE \\u - u h \\ 0 E . 

The bound for the other two differences is obtained by splitting each bilinear form up into its constituent 
terms and applying the definition of polynomial consistency. For the diffusion terms, the polynomial 


15 



consistency property means we consider 


af(ngu, V'i) - a E (n° k u, ipi) = / Kvn^.u • (n^_ x -1) V0i a® 

J E 

= [ (nti-i) (n\/n° k u) • v (v>i - ip) dx+ 

J E 

+ [ (nt, -1) ( K vn°«) • v (v> - n°v>) a® 

J E 

< Ch r+1 \\u\\ r+1E \\u- u h \\ 0 E , 

having applied the Cauchy-Schwarz inequality and the polynomial approximation bounds in the final 
step. The second difference is similarly treated by adding and subtracting terms and applying the 
Caucliy-Schwarz inequality and the polynomial approximation bounds along with the regularity of the 
dual solution ip and the iT-norm error bound 


a E ( u h , H°ip) - a E {u h , U^ip) 



-i) v^-vn^d* 

- I) V(u h - u) • V (ft? — I) ip da;+ 

l_ x -1) vu • v (n? - i)ipdx+ 


+ I V (u h - n° u ) • (ntr -1) (*W) d® 
J E 

< Ch r+1 \\u\\ r+1E \\u - u h \\ 0iE . 


The bounds for the other components of the bilinear form in these differences are treated completely 
analogously. Consequently, we may combine these individual bounds to determine the optimal order 
bound in the statement of the theorem. □ 


6 The Bilinear Forms 

In this section, we introduce a choice of the virtual element bilinear forms a E and b E that satisfy 
the abstract properties presented in Section 3. We remark that the bilinear forms that we pick are 
exactly the same regardless of whether we are considering the conforming or the nonconforming method. 
Moreover, as described in Section 7, the implementation of the two methods differs only in the practical 
construction of the L 2 -projection operators due to the different choice of the degrees of freedom. 

Definition 4. Let E G Th■ A computable (see Definition 1) bilinear form S E : V E /V k (E) x 
V E /V k (E) —» IR is said to be a local admissible stabilising bilinear form if it is symmetric, positive 
definite and it satisfies 

c 0 a E {v h ,v h ) < S E (v h ,v h ) < da E {v h ,v h ) \/v h G V E /V k {E), 

for some constants Cq and C\ independent of E and h. 

Given an admissible stabilising bilinear form S E (-,-) and a computable projection II® : V E — >■ 
Vk(E), we simply define 

a E (u h ,v h ) := (nn^Vu^n^Vv^E + (vn° k u h ,n° k v h ) + s E ((i-n®)u h , (i-nf)v h ), ( 6 . 1 ) 

and 

bh(u h ,v h ) := \ [((3 • n«_ 1 Vt tfc ,Ilgt; h ) - (U° k u h ,f3- II^Vv/O] . (6.2) 

It is clear that any bilinear form A E resulting from (6.1) and (6.2) with S E admissible satisfies Assump¬ 
tion A2, cf. [5]. Moreover, all of the terms in this bilinear form are computable since the construction 
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and degrees of freedom of the space allow us to compute and n°_ 1 Vu^ for any Vh € V E , while S E 
and II® are computable by assumption. The projection II® could be chosen in many ways. A common 
choice would be to use an operator which is already computed, such as the L 2 (.E)-orthogonal projection 
11° or the projection 11)1 used in the definition of the space. Another option could be to choose II® in 
conjunction with the stabilising term in such a way as to incorporate important new features into the 
virtual element method, such as monotonicity, positivity and maximum/minimum principles for the 
numerical solutions, like the mimetic finite difference stabilising terms described in [26]. This topic, 
however, is beyond the scope of the present paper and will be considered for future works. 

We now introduce a choice of admissible stabilising bilinear form. This is essentially the one already 
used in [5]. 

Proposition 1 . Let k e , V • f3 E and 7 E be some constant approximations of k, V • /3 and jj, over E 
respectively. Then, the bilinear form 


1 nE 

S E (v h ,w h ) := ( K E h d E 2 - -V • f3 E h d ~ l + 7 f E h E )^2doi r {vh)doi r (w h ), 

~ r= 1 

for i’h,Wh £ V E /Vk(E), is admissible. 

Proof. The stabilising term S E is an inner product over the finite dimensional space jR/ l ®- dim C p fc(- E )) 
of vectors of degrees of freedom, which is isomorphic to V E /Vk{E). Since a E is an inner product on 
V E and thus also on V E /Vk{E), the existence of the constants Co and C\ of Definition 4 follows from 
the equivalence of the norms induced by these inner products. 

The fact that these constants are independent of h is due to the fact that S E scales the same as 
a E . It is clear that the H 1 part of a E scales like K E h d ~ 2 while the L 2 term ((7 — • (3)vh,Wh) E 

scales like ~p E h d — \ V • j3 E h d ~ 1 . Then, since the degrees of freedom are specifically chosen to scale like 
1 (cf. [7]), the coefficient at the front of S E ensures that the term has the correct scaling even when 
one of the coefficients k, V • /3 or n locally degenerates. □ 

Remark 3. The stabilising term in Proposition 1 is just one of a family of admissible stabilising 
terms, defined as appropriately scaled inner products on the subspace of the degrees of freedom relating 
to functions in V E /Vk{E). Here we have chosen the Euclidean inner product for simplicity. 

6.1 The Effects of Numerical Integration 

In any practical implementation of the method, the coefficients n,/3, and p., have to be approximated, 
meaning that the polynomial consistency properties of Assumption A2 will hold in an approximate way 
in general. One possibility, which we assess here, is to utilise numerical quadratures. Crucially, such 
numerical quadratures will only affect the consistency terms. These have precisely the same structure of 
standard finite elements terms (integral products of polynomial trial and test functions weighted by the 
coefficients), and as such the variational crime introduced by their approximation can be assessed using 
the classical finite element analysis. Indeed, within this section we show that the actual implemented 
methods retain the stability and optimal accuracy properties of the theoretical virtual element methods 
proposed above, provided that the quadrature scheme used is of at least polynomial order 2k — 2. We 
emphasise that this is exactly the same requirement as for the classical finite element methods used to 
solve the same problem (cf. [16]). 

In more concrete terms, suppose that we are approximating integrals over the element E using a 
quadrature rule Qff of degree to, so 


/ gdx « QZ(g) ■= yZ^eg(qe), 
Je t= 1 


for a finite set of quadrature points {qt} E = i and associated weights {w}^ =1 . Thus in practice, the 
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implementation of the method will be based on the perturbed bilinear form 


L 

<*h( u h’ v h) := n(q e )(U° k _iVu h )(q e ) ■ (U° k _ 1 Vv h )(q ( ) 
£=1 


L 

+ “e v(qe)(n 0 kU h )(qe) {U° k v h ){q e ) 
t=\ 


(6.3) 


1 n -E 

+ ( K E h d E - 2 - - + Jl E h d E ) dofr((I -n ®)u h ) dof r ((i -n® )v h ), 

r =1 

and similarly for the skew-symmetric part b/. Once more, we note that the use of quadrature only af¬ 
fects the consistency term. The following two theorems show that the use of an appropriate quadrature 
rule does not affect either the stability or the accuracy of the method. 

Theorem 7. Suppose the quadrature scheme Q ® with m > 2k — 2 has strictly positive weights and is 
exact for the space V 2 k- 2 (E) and/or the set {qt}/ =1 of quadrature points contains a Vk-i (E) unisolvent 
subset. Let 21/, denote the bilinear form Ah with the polynomial consistency integrals approximated using 
the quadrature scheme Qff. Then, 21/, satisfies the stability property of Assumption A2. 

Proof. We first wish to show that the bilinear form a/ defines a norm on V E , or on V///Vo (E) when 
p = 0, equivalent to the norm imposed by a E in either case. Suppose p / 0. Then a/ is clearly already 
a semi-norm and all that remains to be shown is that a E (vh, Vh) = 0 => Vh = 0 . 

Let a E (vh,Vh) = 0. Then, we must have 

QZiin^Vvh) ■ (kHI^X/vh)) =0. 

By the assumptions on Qf/ and the strong ellipticity of k (see (2.2)), this implies that n°_ 1 Vt>/, = 0, 
and consequently we may deduce that either (a), Vh £ Vo(E) or (b), Vh £ V E /Vk{E) where, with a 
slight abuse of notation, we associate V E /Vk{E) with the non-polynomial subspace of V E . 

Suppose case (a) holds, so Vh £ Vo(E). Then we also have 

0 = Qm(M(n°u/,) 2 ) = Q^{pvl), 

since n ° is the identity on Vo(E), and thus we deduce that Vh = 0. 

Alternatively, suppose that case (b) holds, so v k € V E /Vk{E). Then, since a E {vh,Vh) = 0, it 
follows that 

5 B ((i-n®K,(i-nfK) = o, 

and we may deduce that Vh = 0 . 

From this, we may conclude that (a/ (-,-)) 5 is a norm on V E . Moreover, (a B (-,-))5 is also a norm 
on V E and since this is a finite dimensional subspace of iL 1 (£'), the resulting norms are equivalent. 
As with the bilinear form a E , the constants in the equivalence are independent of h due to the correct 
scaling of a E . 

On the other hand, when p = Owe find that (af (•, -)) 2 and ( a E ( •, -)) 5 are both norms on V E /Vo{E) 
and both zero on Vq(E). Again, the fact that V E /Vq(E) is finite dimensional allows us to deduce that 
the two norms are equivalent. 

The stability property for b E also holds because b E (vh, Vh) = b E (vh, Vh) = 0 and it is straightforward 
to check that b E (u h ,v h ) < C||/3|| 00 ||u/ 1 || 1 ,,eIMI 1 ,.e- d 

The next result addresses the questions about the accuracy of the method when the quadrature 
is employed to evaluate the integrals of the virtual bilinear forms, and should be compared with 
Theorem 5. 

Theorem 8. Suppose that Assumptions A1-A3 are satisfied. Let k > 1 be a positive integer and let 
u £ H s+1 (fl) be the true solution to problem (2.4) for some positive integer s. Define r = min(/c, s) 
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and suppose that the coefficients n,(3,p £ W r+1 ’°°(Q), satisfying (2.2) and (2.3). Suppose that the 
right-hand side function f £ is approximated by fh := n^ ax ^ fc _ 2 0 ^/. 

Suppose that the quadrature scheme Q ^ with m > 2k — 2 is exact for the space 'P 2 k- 2 (E). Let 
3ft denote the virtual element bilinear form Ah obtained by approximating the integrals using the 
quadrature scheme Qff and u/, £ Vh be the solution obtained from this scheme. Then, there exists 
a positive constant C, independent of h and u such that 


ll« - «fc|li,h < Ch r ( 1 + h)\\u\\ r+1 + CTil/IU,. 

Proof. Expanding as in the proof of Theorem 2, it may be shown that the only extra term depending 
on the quadrature scheme which arises in the abstract error bound is 

E A%(u n ,w h ) - 3ft {u v ,w h ) 

SU P -ii-ii-> 


where u n := II^m. As usual, we split this term into the different components of the bilinear form 
and bound them separately. We give here the bound for a^{u^,Wh) — <tjf(u,r, w h)', the bound for the 
skew-symmetric term follows analogously. Since the stabilising term is unaffected by the quadrature, 
we only need to bound 


sup 

w h£V* 


Qm(( K Vu-n-) ■ n^jVfflft) - [ (kVUk) ■ n£_i Vw h d® 
_ J_E _ 

1 ,E 


+ sup 



\\ w h || 


Arguing as in Theorem 4.1.4 in [16] and using the stability of the L 2 (E) projector, we find that 


Qm(( K Vu w ) • n l^Vwh) - / {nS7u n ) • n^Vwft da: 

J E 


< Ch r \u\ r+1 , E \\Wh\\i ei 


and 


Qmi^Un) • n ° k w h ) - / {nu^Ul wh da: 


< Ch r \u\ rE \\w h \ 


1 ,E' 


The theorem then follows by treating the skew-symmetric term similarly and combining the result 
with the original iJEnorm bound of Theorem 5. □ 

A similar analysis can be carried out to control the error in the numerical approximation fh of the 
forcing function /, and to recover the optimal order of convergence in the T 2 -norm, see always [16]. 


7 Implementation 

To fix the definition of the space V®, we must first define the computable projection II). : W k / ~ —> 
Vk(E). The original approach in [1] is to contract a particular projector, in that case an elliptic 
projection therein denoted by 11]?. Here instead we first construct the general family of projections 
and subsequently fix a simple choice. 

A projection n) : W)f /A —> Vk(E) may be defined using any computable inner product Be(-, •) on 
Wh / ~Vk{E). For any Vh £ W k /~, we define the polinomial n)uft £ Vk{E) as the solution of 

B E (n* k v h ,m a ) = B E (v h ,m a ) \/m a £V k (E). (7.1) 

Let WJS the Lagrangian basis functions of W ®/~ with respect to the original degrees of 
freedom in Definition 2, and define the matrix D such that Dj„ = dof,(m Q ). Since W k /~ is finite 
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dimensional, B E {‘ v) can be written as the symmetric positive definite matrix B = (B E {fpi,ipj))- From 
this, it can also be seen that (S_E(m a , V’i)) = D 1 B and (BE(m a ,mf 3 )) = D 2 BD. 

Define the action of II£ on the shape functions {?/>*}"=! through the matrix II®, where 

Nd,k 

If-fcV’i = ^ ' ^a(Hfc )ait 
oc—1 

so that the j-th column of II® contains the coefficients of the expansion of the polynomial Tl k ipj in the 
monomial basis {m a }^=i- Consequently, the projection problem (7.1) can be written in the matrix 
form 

nf = (d t bd) _1 d t b 

The matrix D 2 BD is invertible due to the assumption that B E (•, •) is an inner product on W)f/~- 
Finally, we have to choose the bilinear form B to fix the projection. This is equivalent to picking 
any symmetric positive definite matrix B. Here we choose B = I, yielding the simple choice 

nf = (D t D)” 1 D 2 '. (7.2) 

Having thus chosen the projection n£ and therefore the space V k , we can introduce the Lagrangian 
basis j of V k associated with the degrees of freedom in Definition 2. These shape functions are 

necessarily different for the conforming and the nonconforming methods, although the implementation 
of the two methods is formally the same since we are only concerned with the degrees of freedom. 

We shall determine the local matrix A k (fa, c/)j) and right-hand side vector (fh^^E associated to 
the Lagrangian basis introduced above. To this end, we first need to evaluate the projections 
and for all * = 1 ,..., n E - 

The polynomial n^,; is the solution of the projection problem: 

K,H° t (ii) J? = K<),)£ Va = l,...,N d ' k . (7.3) 

Since H k fa is both a polynomial of degree k and a function in the virtual element space V k , it can be 
expanded on the monomials generating Vk(E) and the shape functions generating Vj® as 

Nd,k TIE 

H 14* = E m «(n° k ) ai = 

a=l j=l 

and the coefficients of these expansions are collected in the matrices 11° and 11°’^, respectively. The 

matrix Iljj*’^ will be used at the end of the subsection to compute the stabilising term. By comparison, 
it follows that 11°’^ = DII°. 

We reformulate the projection problem (7.3) in matrix form as HII° = C, where the coefficients of 
matrices H and C are given by 

Cm = (m a ,<j>i) E and H a/3 = (m a mp) E 

for a = 1,..., Nd,k and i = 1,... ,n E . Since the space V)® does not use the extra degrees of freedom in 
Definition 3, the matrix C must be constructed in two parts according to the definition of the space. 
Then, we have 


r — 


(m a ,<j>i) E if m a G M k - 2 {E), 
(m a ,A* k fa) E if m a G Ml_faE) U M* k (E). 


(7.4) 


With the choice of nji presented above, C becomes 


r . — 

'-m, — 


{m a ,<pi) E if rn a G M k - 2 (E), 

(H(D T D) -1 D T ) ai if m a G M k _ x {E) U M k (E), 


(7.5) 
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where the moments of fa are simply degrees of freedom, so C is fully computable. 

The other crucial term which must be computed is \7 fa, where the projection is defined 
componentwise such that 

(II l-iVfa,m a )E = (V fa,m a ) E Vm„ g ( M k -i(E )) 2 

= / m a ■ nfads - (fa, V • m a ) E . 

JdE 

The second term on the right-hand side of this expression is simply a combination of the internal 
degrees of freedom of fa. For the nonconforming method, or for the conforming method when d = 2, 
the first term is also just a combination of edge degrees of freedom of fa. However, for the conforming 
method when d = 3, we must compute this term using the L 2 -ortlrogonal projection of fa on each face 
s C dE as 

/ m a ■ nfads = ^ / m a ■ nll^fads. 

J dE sCdE^ s 

The projection n°’ s can be computed on each face of E exactly as when d = 2. 

Thus, we end up with d linear systems to compute the projection of the d components of fa, 

namely = R xe where G a p = {m a ,mp)E with m a ,mp € M. k -i(E) and 

( R xe )ai = Y] / m a n e fa ds - (fa, ——)e, 

sUe^ 3X l 

for l = 1,..., d, where ng is the f'-th component of n. Then, we have that 

Nd,k -1 


H 


0 

k ~ x dxi 


= ^o( n fc-l) 


CK— 1 


Having the matrices 11° and n l = l,...,d, we are able to easily implement all the terms in 
the local bilinear form AF. Indeed, for the term of aF that contains the diffusion coefficient K\, we 


have: 


(“ 


tto dfa o 
Un ^dxi k ^ 1 dx n >E 


Nd,k-1 

= ( K Ur. m *, m i3) E ( n k-[) m ( U k-l) 

a,/ 3=1 


Pj 


, n = 1 ,..., d. 


For the reaction term in af, we easily find that: 


Nd,k 


(niL 0 k fa,n 0 k <t>j) E = (/ im «, m ^(n“) (1 ,(n2) w . 


a,/3—l 


For the skew-symmetric bilinear form bf, first notice that / 3 • n°_ 1 V</>j = Ylt=i AH°_ 1 |^,,:. where fa 
is the Z-th component of (3. Therefore, the first term of bf is given by: 


IE 


d Nd,k Nd,k -1 

(3U° k fa • H^V^- d* = £ £ E (fa m a , mp) E (nl) a .(n°£\) 

1=1 a-1 0—1 




and a similar expression is found for the second term by exchanging i and j. 

Finally, according to the expressions given in Proposition 1, the stabilising term is given by 

Sf E ((i-ng) 0 i, (i-n°)^) = ( k E h d E - 2 - \Vfa3 E h d ~' +j E h d E ) ((i - n°/) T (i - n°-*)).. 

since dof r ((I-H°)^) = (l-n^). r . 
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Similarly, since we have the projector 11°. at our disposal, we might as well use fh := II°/|b to 
approximate f\ e- In this case the right-hand side vector is given by: 


{fhAj)E = 


IE 


n kffa dx 



(7.6) 


In practice, the L 2 -products on the right-hand side of the above formulas have to be somehow approx¬ 
imated in accordance with the theory presented in Section 6.1. 


8 Numerical Results 

All the numerical experiments presented in this section are obtained using (7.6) for the approximation 
of the right-hand side, the choice given by (7.2) for n£, setting Ilf = 11°, and using (7.5) for the 
definition of matrix C. However, a comparison with the implementation using nf = nf and = nf 
(cf. [1]) in (7.4) did not reveal any significant difference in the behaviour of the method. 

The numerical experiments are aimed to confirm the a priori analysis developed in the previous 
sections. In a preliminary stage, the consistency of both the conforming and nonconforming VEM, i.e. 
the exactness of these methods for polynomial solutions, has been tested numerically by solving the 
elliptic equation with boundary and source data determined by u(x , y) = x m + y m on different set of 
polygonal meshes and for to = 1 to 4. In all the cases, we measure an error whose magnitude is of the 
order of the arithmetic precision, thus confirming this property. 

To study the accuracy of the method we solve the convection-reaction-diffusion equation on the 
domain f 1 =]0, l[x]0,1[. The variable coefficients of the equation are given by 

, , ( 1 + y 2 — xysin(27ra;) sin(27ru) 

v y y —a:ysm(27ra;) sm(27ry) 1 + x 2 

ft x > y) = ( 3 ( 3 tf-ly + 3) )’ 7 ( *’ y) = x2 + y 3 + 1 - 

The forcing term and the Dirichlet boundary condition are set in accordance with the exact solution 

u(x, y) = sin(27ra) sin(27n/) + x 5 + y 5 . 

The performance of the methods presented above are investigated by evaluating the rate of conver¬ 
gence on three different sequences of five meshes, labeled by M\, M 2 and M 3 , respectively. The top 
panels of Fig. 5 show the first mesh of each sequence and the bottom panels show the mesh of the first 
refinement. The meshes in A4i are built by partitioning the domain Q into square cells and relocating 
each interior node to a random position inside a square box centered at that node. The sides of this 
square box are aligned with the coordinate axis and their length is equal to 0.8 times the minimum 
distance between two adjacent nodes of the initial square mesh. The meshes in A 42 are built as follows. 
First, we determine a primal mesh by remapping the position ( 2 :, y) of the nodes of an uniform square 
partition of fi by the smooth coordinate transformation [24]: 

x = x + (1/10) sin(27rir) sin(27ry), 
y = y + (1/10) sin(27ra;) sin(27ry). 

The corresponding mesh of M 2 is built from the primal mesh by splitting each quadrilateral cell into 
two triangles and connecting the barycenters of adjacent triangular cells by a straight segment. The 
mesh construction is completed at the boundary by connecting the barycenters of the triangular cells 
close to the boundary to the midpoints of the boundary edges and these latters to the boundary vertices 
of the primal mesh. The meshes in .AI 3 are obtained by filling the unit square with a suitably scaled 
non-convex octagonal reference cell. 

All the meshes are parametrised by the number of partitions in each direction. The starting mesh 
of every sequence is built from a 5 x 5 regular grid, and the refined meshes are obtained by doubling 
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Figure 5: First (top) and second (bottom) mesh of the three mesh families. 


this resolution. Mesh data for each refinement level, i.e., numbers of mesh elements, number of edges, 
number of vertices, are reported in Table 1. 

Approximation errors are measured by comparing the polynomial quantities 11° Uh and , Viq,, 
which are obtained by a post-processing of the numerical solution, with the exact solution u and 
solution’s gradient V«. The relative errors for the approximation of solution u and its gradient in 
function of the mesh size h are shown in the log-log plots of Fig. 6 for the mesh sequence M\, Fig. 7 
for the mesh sequence M 2 , and Fig. 8 for the mesh sequence M 3 . The values of the measured error 
are labeled by a circle for the conforming VEM and by a square for the nonconforming VEM. The 
plots on the left show the relative errors for the approximation of the solution, while the plots on the 
right show the relative errors for the approximation of the solution’s gradient. The expected slopes are 
shown for each error curve directly on the plots. The numerical results confirm the theoretical rate of 
convergence. The conforming and nonconforming VEMs provide very close results on any fixed mesh, 
with the conforming method slightly over performing the nonconforming VEM in few cases. Similar 
results (not shown) are observed when comparing the two methods with respect to the respective 



Randomised quadrilaterals 

Remapped hexagons 

Non-convex octagons 

n 

Me 

Ms 

M v 

h 

Me 

M s 

M v 

h 

Me 

Ms 

M v 

h 

1 

25 

60 

36 

0.331 

36 

125 

90 

0.328 

25 

120 

96 

0.291 

2 

100 

220 

121 

0.186 

121 

400 

280 

0.185 

100 

440 

341 

0.146 

3 

400 

840 

441 

0.094 

441 

1400 

960 

0.097 

400 

1680 

1281 

0.073 

4 

1600 

3280 

1681 

0.047 

1681 

5200 

3520 

0.049 

1600 

6560 

4961 

0.036 

5 

6400 

12960 

6561 

0.024 

6561 

20000 

13440 

0.025 

6400 

25920 

19521 

0.018 


Table 1: Mesh data for the meshes in Mi, M 2 , and M 3 ; Me, M s and M v are the numbers of mesh 
elements, interfaces and vertices, respectively, and h is the mesh size parameter. 
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Figure 6: Error curves for the conforming VEM (circles) and the non-conforming (triangles) applied to 
the mesh family of randomised quadrilaterals with k = 1, 2, 3, and 4. The left panels show the relative 
L 2 error; the right panels show the relative H 1 errors. The expected slopes are indicated by triangles. 


number of degrees of freedom. Indeed, for each mesh shown, the difference on the number of degrees 
of freedom does not depend on the polynomial degree k and is about equal to the number of elements, 
in favour of the conforming VEM. 


9 Conclusion 

We have introduced a unified abstract framework for the Virtual Element Method, through which 
conforming and nonconforming VEMs for solving general second order elliptic convection-reaction- 
diffusion problems with non-constant coefficients in two and three dimensions are defined, analysed, and 
implemented in a largely identical manner. We have shown that both methods produce solutions which 
converge to the true solution at the optimal rate in the H 1 - and L 2 -norms, supported by numerical 
experiments on a variety of different mesh topologies including non-convex polygonal elements. 

The framework is based on assuming that the L 2 -projector onto the polynomial subspace of the 
virtual element space is computable, in this respect following the approach of [1, 8]. By generalising the 
process considered in [1], we have introduced families of new possible conforming and nonconforming 
virtual element spaces in which the L 2 -projection is indeed exactly computable directly from the 
degrees of freedom used to describe the space. From this family we have detailed a particular space 
for which the implementation of the L 2 -projection takes a simple form independent of the method, 
the polynomial degree, and the space dimension. It also becomes apparent that since the accurate 
approximation of the problem’s data is only needed to evaluate the polynomial consistency part of the 
bilinear form, the variational crime theory classical of finite element methods applies to the virtual 
element setting. Extensions of the present framework to include stabilisation techniques for convection- 
dominated diffusion problems and the design of virtual element methods for Stokes problems will be 
considered in future works. 
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